Loading libraries

library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
library(scales)
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(leaflet)

Read in data

ak_regions <- read_sf("data/ak_regions_simp.shp")

#ploting with base R
plot(ak_regions)# something is happening realted to the coordinate reference system

#what is ak_regions?
class(ak_regions)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
#it is a dataframe so we can call head()
head(ak_regions)
## Simple feature collection with 6 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.2296 ymin: 51.15702 xmax: 179.8567 ymax: 71.43957
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 4
##   region_id region     mgmt_area                                   geometry
##       <int> <chr>          <dbl>                         <MULTIPOLYGON [°]>
## 1         1 Aleutian …         3 (((-171.1345 52.44974, -171.1686 52.41744…
## 2         2 Arctic             4 (((-139.9552 68.70597, -139.9893 68.70516…
## 3         3 Bristol B…         3 (((-159.8745 58.62778, -159.8654 58.61376…
## 4         4 Chignik            3 (((-155.8282 55.84638, -155.8049 55.86557…
## 5         5 Copper Ri…         2 (((-143.8874 59.93931, -143.9165 59.94034…
## 6         6 Kodiak             3 (((-151.9997 58.83077, -152.0358 58.82714…
#to see the coordinate system
st_crs(ak_regions)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Transforming coordinate systems

ak_regions_3338 <- ak_regions %>% 
  st_transform(crs = 3338)

plot(ak_regions_3338)

We can use tidyverse functions to the sf object

ak_regions_3338 %>%
  filter(region == "Southeast") %>% #select rows
  select(region) # select column
## Simple feature collection with 1 feature and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 559475.7 ymin: 722450 xmax: 1579226 ymax: 1410576
## epsg (SRID):    3338
## proj4string:    +proj=aea +lat_0=50 +lon_0=-154 +lat_1=55 +lat_2=65 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
## # A tibble: 1 x 2
##   region                                                           geometry
## * <chr>                                                  <MULTIPOLYGON [m]>
## 1 Southeast (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7, 12967…

EPSG CODES

epsg.io interesting link

3338- Alaska Albers 4326- WGS84 (GPS) 3847- pseudo mercator (Google maps, Open street maps)

Spatial joins and summaries

pop <- read.csv("data/alaska_population.csv", stringsAsFactors = F)

head(pop)
##   year     city      lat       lng population
## 1 2015     Adak 51.88000 -176.6581        122
## 2 2015   Akhiok 56.94556 -154.1703         84
## 3 2015 Akiachak 60.90944 -161.4314        562
## 4 2015    Akiak 60.91222 -161.2139        399
## 5 2015   Akutan 54.13556 -165.7731        899
## 6 2015 Alakanuk 62.68889 -164.6153        777
class(pop)#it is a data frame and we need a sf object
## [1] "data.frame"
#we are going to convert a data frame to a sf

pop_4326 <- st_as_sf(pop,
                     coords = c("lng", "lat"), # first is x and then y
                     crs = 4326,
                     remove = F) #It keeps the lat long coordinates we originally have

Making the join

#pop_joined <- st_join(pop_4326, ak_regions_3338, join = st_within)

we get an error because the projections are different

We transform our pop

pop_3338 <- pop_4326 %>% 
  st_transform(crs = 3338)

We join

pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)

We plot it

plot(pop_joined)

we are going to calculate population by region

pop_region <- pop_joined %>% 
  group_by(region) %>% 
  summarise(total_pop = sum(population))

droping geometry column. For help: ?sf::tidyverse

pop_region <- pop_joined %>% 
  as.data.frame() %>% 
  group_by(region) %>% 
  summarise(total_pop = sum(population))

now we join the geometry from ak_regions to pop_region

pop_region_3338 <- left_join(ak_regions_3338, pop_region, by = "region")

plot(pop_region_3338)

calculation total population on management area (mgmt column)

pop_mgmt_3338 <- pop_region_3338 %>% 
  group_by(mgmt_area) %>% 
  summarise(total_pop = sum(total_pop), do_union = FALSE) #do.union= F keeps the lines of the regions

plot(pop_mgmt_3338["total_pop"])

write_sf(pop_region_3338, "data/ak_regions_pop.shp", delete_layer = TRUE)# delete layer it removes previous

Make maps!

using scale package to change labels to comma

ggplot() +
  geom_sf(data = pop_region_3338, aes(fill = total_pop)) + 
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma) #labels comes from scale packages.

adding points

ggplot() +
  geom_sf(data = pop_region_3338, aes(fill = total_pop)) + 
  geom_sf(data = pop_3338, aes(), size = 0.5) +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma) #labels comes from scale packages.

adding rivers

rivers_3338 <- read_sf("data/ak_rivers_simp.shp")

st_crs(rivers_3338)
## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=aea +lat_0=50 +lon_0=-154 +lat_1=55 +lat_2=65 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"
ggplot() +
  geom_sf(data = pop_region_3338, aes(fill = total_pop)) + 
  geom_sf(data = rivers_3338, aes(size = StrOrder), color = "black") +
  geom_sf(data = pop_3338, aes(), size = 0.5) +
  scale_size(range = c(0.01, 0.2), guide = F) + #to avoid problems with different scales
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma) #labels comes from scale packages.

Getting basemaps using ggmap it works with 3857

pop_3857 <- pop_3338 %>% 
  st_transform(crs = 3857)

We are getting a stamenmap from stamenmap

# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))
  
  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
  
  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}
bbox <- c(-170, 52, -130, 64)

ak_map <- get_stamenmap(bbox, zoom = 4)
## Source : http://tile.stamen.com/terrain/4/0/4.png
## Source : http://tile.stamen.com/terrain/4/1/4.png
## Source : http://tile.stamen.com/terrain/4/2/4.png
## Source : http://tile.stamen.com/terrain/4/0/5.png
## Source : http://tile.stamen.com/terrain/4/1/5.png
## Source : http://tile.stamen.com/terrain/4/2/5.png
ak_map_3857 <- ggmap_bbox_to_3857(ak_map)

class(ak_map_3857)
## [1] "ggmap"  "raster"

Mapping

ggmap(ak_map_3857) + 
  geom_sf(data = pop_3857, aes(color = population), inherit.aes = F) +
  scale_color_continuous(low = "khaki", high =  "firebrick", labels = comma)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Interactive mapping using leaflet

Here we define a leaflet projection for Alaska Albers, and save it as a variable to use later.

epsg3338 <- leaflet::leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:3338",
  proj4def =  "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
  resolutions = 2^(16:7))

let’s use st_transform yet again to get back to WGS84

pop_region_4326 <- pop_region_3338 %>% 
  st_transform(crs = 4326)

generate the map

pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1,
                    label = ~region) %>% 
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")
m

with ti

m <- leaflet() %>%
  addTiles() %>% 
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1,
                    label = ~region) %>% 
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")
m